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Abstract 



The ever-expanding catalog of detected super-Earths calls for theoretical stud- 
ies of their properties in the case of a substantial water layer. This work considers 
such water planets with a range of masses and water mass fractions (2 — 5 M Earth , 
0.02 — 50% H 2 0). First, we model the thermal and dynamical structure of the 
near-surface for icy and oceanic surfaces, finding separate regimes where the 
planet is expected to maintain a subsurface liquid ocean and where it is expected 
to exhibit ice tectonics. Newly discovered exoplanets may be placed into one 
of these regimes given estimates of surface temperature, heat flux, and gravity. 
Second, we construct a parametrized convection model for the underlying ice 
mantle of higher ice phases, finding that materials released from the silicate - 
iron core should traverse the ice mantle on the time scale of 0.1 to a hundred 
megayears. We present the dependence of the overturn times of the ice mantle 
and the planetary radius on total mass and water mass fraction. Finally, we 
discuss the implications of these internal processes on atmospheric observables. 

Subject headings: convection, planetary systems 
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Introduction 



The discoveries of planetary systems around other stars has now led us into the realm 
of planets in the mass range from about 2 to 10-15 Me ar ^, which are expected to be 
terrestrial in nature and unlike the gas giants in our Solar System. This new type of planets 



has been collectively called "super-Earths" (Melnick et al. 2001). Initially, mass is the 



main physical parameter available. As such, observers who have discovered nearly a dozen 



such planets to date (e.g., Mayor et al. 2009 and references therein) cannot yet distinguish 



between the various interior structures that are expected in this mass range, i.e. water-rich, 
dry, or gas-dominated ('mini Neptunes') types. This is expected to change shortly, with 



the first super-Earth discovered in transit, CoRoT-7b (Leger et al. 2009), and many more 
anticipated to come from NASA's Kepler mission, affording derivation of mean densities 
and bulk composition. 

Theoretical models of planet formation anticipate large numbers of water-rich 



super-Earths (Ida and Lin 2004, Raymond et al. 2004, Kennedy et al. 2006). In this 



paper we define as water planet any planet in the super-Earth mass range composed of 
> 10% H 2 by mass around a silicate - metal core and that lacks a significant gas (H and 
He) layer. 



Water planets were proposed by Kuchner (2003) and Leger et al. (2004), and modeled 



with detailed T-dependent equations of state (EOSs) by Valencia et al. (2007a), Sotin 



et al. (2007), and Grasset et al. (2009). Mass-radius relations have been computed to 



help distinguish between the different types of super-Earths when observations become 



Valencia et al. 


2007b, 


Fortney et al. 


2007, 


S eager et al. 


2007 



Selsis et al. 2007). Ehrenreich and Cassan (2007) and Tajika (2008) considered cold, 



ice-covered water planets. With the prospect of high-quality observational data on water 
planets becoming available (CoRoT, Kepler, JWST), there is strong motivation to take the 
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modeling described above to a new level and focus on the geophysical details of the water 
layer (s), which is the goal of this work. 

This work begins by constructing hydrostatic models that take into account the effect 
of ice phase transitions up to Ice X. We model nine ocean planets with total mass ranging 
from 2 to 10 M Earth and water mass percents from 0.02% (Earth-like) to 50%. Section [2] 
explains the creation of these models. 

Section [3] presents a pair of models to constrain the thermal and dynamical structure 
within the H2O layer of these planets. We consider a variety of imposed conditions such as 
surface temperature (above and below freezing) and heat flux. The first model governs the 
behavior of the near-surface region and considers both conduction and solid-state convection 
as modes of heat transport. The second is a parameterized model for the underlying mantle 
of higher-phase ices. The results from these models are presented in Section |4j 

We follow with a discussion (Section [5j of remaining uncertainties and of the observable 
implications of these thermal models, especially regarding the transport of material from 
the silicate - metal core to the surface. 



2. Hydrostatic Models 

We construct hydrostatic models of nine hypothetical planets with total masses of 2, 
5, and 10 M.Earth- For each mass value, we consider a dry, Earth-like (simplified to be 32% 
Fe, 68% MgSi02i and 0.02% free H 2 0) case and two cases with 25% and 50% water by 
mass around a core with Earth-like composition. These planets are simplified to contain up 
to five distinct compositional layers: iron, MgSiO^ perovskite, Ice X, Ice VII, and Ice Ih. 

We use the canonical equations for a spherical body in hydrostatic equilibrium in 
conjunction with an EOS for each material to obtain the preliminary dimensions of the 
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exoplanets. We follow the results of Seager et al. (2007) and choose the Vinet and 



fourth-order Birch-Murnaghan EOS (Vinet and Ferrante 1987, Shankar et al. 1999) for 



iron and MgSiO^ perovskite, respectively. We choose the third-order Birch-Murnaghan 
EOS for all three water ice phases. Our choices of EOS and parameters used are summarized 
in Table [T] The effect of the choice of EOS is minor for our purpose here (e.g., Valencia 
et al. |[2007a| Seager et al. 2007, Grasset et al. 2009). Following these studies, we use 
isothermal EOSs for all regions. 

Several approximations are employed here. Clearly, the pure chemical composition 
assumed in each of the three regions is non-physical, but its associated error is acceptable. 
For instance, the replacement of 25% of the Mg by Fe in the post-perovskite MgSiOs 



structure results in a 2.5% increase in density (Caracas et al. 2008). Any changes in the 



property of the silicate mantle due to different internal water content compared to that of 
the Earth are neglected. 

As for phase changes, we neglect any transitions in the Fe core. Due to the lack of 
necessary parameters for higher phases, the entirety of the silicate mantle is assumed to 
be in the perovskite phase of MgSiO^. The higher-pressure post-perovskite phase is stable 



from around 125 GPa up to a currently unknown phase transition (e.g. Mashimo et al. 



2006) or to its dissociation into MgO and S1O2 at about 1000 GPa (Unemoto et al 2008) 



making it the predominant component of the silicate layer of all of our modeled planets 
except the 2 M^arth, Earth-like one. Even so, we find that this as well as the above described 
approximations are acceptable, given that our focus in this work is on the overlying water 
layers. 



Following the experimental phase diagram for high pressure ice phases (Song et al. 



2003, Lobban et al. 1998, Goncharov et al. 1999), we use the elastic parameters of Ice VII 



between 2.2 and 60 GPa, while Ice X is assumed to exist above 60 GPa. In fact, no static 
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laboratory experiments constrain the behavior of ice at pressures above 170 GPa, while the 
pressure at the silicate-ice boundary in the 5 and 10 MEarthi 50% water planets reaches 



212 and 437 GPa, respectively. Recent experiments by Lobban et al. (1998) have led to 
the possibility of a higher ice phase at 150 GPa and above. However, for lack of further 
experimental data, we approximate all potential higher phase ices as Ice X. 

In our hydrostatic models, we ignore the possible existence (depending on thermal 
profile) of other phases of water (e.g. liquid, Ice III, Ice VI) between the Ice Ih and Ice VII 
regions. Due to the relatively small thickness of these near-surface water layers, the effect 
of these approximations on key parameters such as surface gravity g s and planetary radius 
r p is negligible. 

A scaled diagram of the H 2 region of our 5ME ar th, 50% H 2 planet is shown in 
Figure [1} Note that the details in the Ice VI, liquid water, and Ice Ih layers were found 
after the thermal modeling in the subsequent sections. 



3. Thermal Models 

We include two sets of thermal models. The first models the near-surface region of 
the planet, both in the case of a frozen surface and a liquid surface. The second models 
convection in the underlying mantle of higher-phase ices. 



3.1. Near surface convection 



We consider the case where the surface temperature is below the freezing point. We 
model both conductive and convective heat flow in the near-surface region. We define the 
superficial ice layer where heat transport is solely due to conduction as the "crust." For the 
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thickness of the crust, we use the simple expression for a conductive thermal profile: 



D 



k(T B - T s 

qs 



X 



Where D is the thickness of the crust, Tg and Ts are the temperatures at the base 
and at the surface, qs is the surface heat flux, and k is the thermal conductivity of Ice Ih, 



equal to 3.3 Wm 1 K 1 (Spohn and Schubert 2003). Ts and qs are left as the independent 



variables in this analysis. The former is determined by the planet's insolation, albedo, 
and atmospheric effects, while heat flux can be estimated by scaling the planetary heat 
output with the total mass of silicates and metals (M e ff) compared to the Earth, and then 
adjusting this quantity to the desired radius (r): 



M, 



Qs — QEarth- 



eff 



(2) 



Where q B 



arth 



^Earth r Earth 

0.087 W/m 2 ( |Turcotte and Schubert l2002[ ). This estimate for heat flux 



is only used in this work for analyzing the convection of the ice mantle (Section 3.2). 



We evaluate the convective stability of the surface ice layer by finding its Rayleigh 
number (Ra) and comparing it to a critical Rayleigh number (Ra crit ). The former is given 
by: 



Ra 



gapD 3 (T B -T s 



(3) 



ktj(Tb,T s ) 

Where g is the surface gravity, a = 1.6 x 10 -4 K~ x is the thermal expansivity of Ice Ih, 



p = 917 kg m 3 is its density, and k = 3.7 x 10 6 m 2 s 1 is its thermal diffusivity (Spohn 



and Schubert 2003). The viscosity 77 is strongly temperature dependent. Given the low 



pressure regime found in the crust, we assume a Newtonian, volume diffusion-dominated 
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rheology following the work of Mueller and McKinnon (1988): 



77 (T) = ?7oexp[ 



Q* 

RT„ 



(4) 



We adopt values of 139 Pa-s and 56.7 kJ/mol for 7/0 and Q*, respectively, corresponding 



to a grain size of m 0.2mm for our temperature range (McKinnon 2006). 



We note that the stability of the conductive boundary layer has been used to derive 
the relation Nu = A ■ Ra@, where Nu is the Nusselt number and A and are constants. 
Scaling relations of this form are used to relate the convective heat flux to the parameters 



governing convection (Howard 1966), and they yield the same result as more complete 



boundary layer models (O'Connell and Hager 1980, Turcotte and Schubert 2002, Valencia 



and O'Connell 2009). Where applicable (i.e. in the case of convection in an Ice Ih shell 
overlying a liquid ocean), this scaling has been verified to corroborate our results. 

We evaluate Equations [TJ [3] and [4] under two distinct formulations. For relatively 
high Ts (> 190 — 220 K, depending on curve being calculated), we find that the viscosity 
contrast is relatively small (Vmax/Vmin i$ 10 3 ). We therefore employ the small viscosity 
contrast (SVC) presciption, with characteristic viscosity taken to be the value corresponding 
to the average temperature: T v ^svc — {Tb + Tg)/2. This approximation has been found to 
be appropriate for viscosity contrasts of up to f] m ax/Vmin ~ 10 4 (Dumoulin et al. 1999). 



In this regime, we use a critical Rayleigh number of Ra cr i t ,svc — 2000, corresponding to 
viscosity contrasts of 10 3 (Schubert et al. 2001). 



For relatively low T$ and correspondingly high viscosity contrasts, we use the 
stagnant-lid formulation to describe convection. In this case, convection is limited to 
an adiabatic zone in the lowest part of the ice layer. The characteristic viscosity of the 
convection cell is found by evaluating Equation [4] with the temperature at the bottom of 
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the conductive crust: T V) sl — Tb (Solomatov 1995). The corresponding critical Rayleigh 



number is dependent on the magnitude of the viscosity contrast: Ra cr i tt sL — 20. 9p A where 
p reflects the logarithm of the viscosity contrast between the top and bottom of the full ice 



layer (Solomatov and Moresi 2000): 



Q',T B -T : 



P-^V 1 ) (5) 



Where T is the temperature in the actively convecting region. It is cooler than Tg by 
a small amount on the order of 10 K: (Tb — ~ RTf /Q*). 

Using the appropriate formulation for the given viscosity contrast, we solve for the 
maximum thickness of the conductive crust and the bottom temperature Tb as functions of 
the imposed parameters q$, Tg, and g. If this maximum conductive thickness is greater than 
that of the Ice Ih layer before it meets the melting curve, then our crust directly overlies 
a liquid ocean, which has an adiabatic temperature profile. Otherwise, we find solid-state 
convection immediately below the crust. 

In the case of a liquid surface, we model the temperature gradient in the surface ocean 
to be adiabatic. The thickness of this ocean is found by identifying the intersection of this 
adiabat with the freezing curve of water at depth. 



3.2. The ice mantle 

Regardless of the near-surface structure, convection is always present in the underlying 
layers of higher phase ices (henceforth called the "ice mantle"). We assume for now and 
show below that convection within the ice mantle is not partitioned by phase transitions. 
We also assume a rectangular geometry for the convective layer, a valid approximation 
given that the analysis focuses on only the boundary layers. For the case where the ice 
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mantle underlies a liquid ocean, we adopt a thermal profile as shown in Figure [2] This 
thermal profile is partitioned into an adiabatic region between two boundary layers. The 
Bottom Boundary Layer (BBL) is a conductive zone that corresponds to the thin region 
where the cold material recently arrived from above is slowly heated by conduction from 
the underlying silicates until it rises again. 

On the other hand, the temperature profile in the top boundary layer (TBL) is 
constrained to follow the melt curve. The dynamical characteristic of this region is outlined 
further in the discussion (Section 5). The top interface temperature and depth Zt 
correspond to the temperature and location of the bottom of the liquid ocean and are 
constrained by the above considerations of the near-surface thermal profile. To constrain 
the quantities T , AT a ^, and ATbbl ( see Figure [2]), we construct a convection model that 
uses three independent relations: one governing each boundary layer and one describing 
the adiabat. We note that we do not explicitly use any existing Nusselt - Rayleigh number 
scaling laws. However, our model considers precisely the same physics that forms the 
basis of such scaling laws- namely by assuming that both boundary layers are at the verge 



of convective stability (see Valencia and O'Connell 2009 the Appendix for a detailed 
derivation). 

First, we examine the convective stability of the TBL. In all six of our examined ice 
mantles (of the 25% and 50% H2O planets), when approximating heat flux by Equation [2j 
the TBL is composed entirely of Ice VI. The thickness (Stbl) and bottom temperature (T ) 
in the TBL is therefore related simply by following the melt curve of Ice VI: 

T - Tyj 

Otbl = z T (6) 

m VI 

Where Tyi and myi are the T-intercept and slope of the Ice VI melt curve, respectively, 
and Zt is the depth of the bottom of the liquid ocean. The Rayleigh number of the TBL 
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can then be calculated and set equal to the critical Rayleigh Number (Ra crit ), chosen to be 



1000 to reflect the mild viscosity contrast within the layer (Turcotte and Schubert 2002): 



p - Rn 9ap ( T ° — Tt "> 5 TBL m 

J^O'TBL — ti&crit — { < ) 

KVTBL 

Where f]TBL is the viscosity in TBL evaluated for the average temperature and pressure, 



as appropriate for smaller viscosity contrast (McKinnon 1998). T may be obtained from 



Equations [6] and [7j allowing for the evaluation of the adiabatic temperature gradient via 



the following relationship (Turcotte and Schubert 2002): 



dT{z) ^ T(z)ga(P(z),T(z)) 
dz ~ C p (T(z)) 1 ' 

Where Cp(T) is the isobaric specific heat capacity. Equation [8] can be integrated 
numerically from bottom of the TBL to the top of the BBL to find the temperature change 
in the adiabatic regime AT a ^ as a function of T . 

Finally, we evaluate the Rayleigh number of the BBL and set it to the critical Rayleigh 
number Ra crit , again chosen to be 1000. Defining 5bbl to be the thickness of the BBL, we 
can write the Rayleigh number and the heat flux as: 

R&bbl = Ro-crit = BBL BBL (9) 

KT)BBL 

kAT BBL 

Obbl = (10) 

Qbbl 

Where tjebl is the BBL viscosity, which is dependent on the value of To and T a d, and 
is once again calculated for small internal viscosity contrast. Combining Equations M and 



10 allows the calculation of the final unknown, ATbbl- 
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Due to the much higher pressure and temperatures in the ice mantle, the assumptions 
of Newtonian rheology and constant a as in our near-surface analysis are insufficient. Under 
the stress regime in our ice mantles (several MPa, see below), dislocation creep, which 



has a significant stress dependence, is expected to dominate (Durham and Stern 2001). 



Furthermore, the viscosity change between known ice phases under this stress regime is 
not great (< 10 2 ). Hence, we use a dislocation creep model for the viscosity of Ice VI, the 



highest pressure phase measured to date, for the ice mantle (Durham et al. 1997): 



r)(P, T) = 6.65 x 10 19 o-" 3 - 5 exp[-(£* + PV*)/RT] 



'IV, 



Where a is a characteristic shear stress, R is the ideal gas constant, and E* and V* 
are the activation energy and volume, measured to be 110 kJ/mol and 1.1 x 10~ 5 m 3 /mol, 
respectively. Owing to the lower viscosity of ice compared to silicates, the shear stress is 
approximated to be a value somewhat lower than that of the Earth's mantle, 2 x 10 6 Pa. 
This value is self-consistent with our results. The activation volume V* itself changes 
with depth; it can be approximated to shrink with pressure as a vacancy in the material. 



Following the analysis of O'Connell (1977), the effective bulk modulus of this vacancy is: 



K* = V* 



dP 



2(1 - 2u) 



3(1 



-K 



v) 



(12) 



Where K is the bulk modulus of the material and v = 0.35 is the Poisson ratio. By 



expanding K to first order in pressure and integrating Equation [12j we obtain the activation 
volume as a function of pressure: 



fC* _i_ x'* P i 
V*(P) = VZ( 0+ 



K* +K'*P 



(13) 



Combining Equations 11 and 13 gives viscosity as a function of depth for a specified 
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thermal profile. Under this formulation for viscosity, our results show that maximum 
viscosity contrasts in the TBL are on the order of 10 2 , while that of the BBL are similar 
in most cases, except for the 2 MEarth planets, where the contrasts are between 10 5 and 
10 7 . Despite this, we maintain the use of the small viscosity contrast formulation for BBL 
viscosity due to (1) the fact that varying tjbbl through its entire possible range does not 
change our ultimate results (overturn timescale, etc) by more than a factor of a few and (2) 
the lack of existing theory for bottom boundary layers with high viscosity contrast. Note 
that the viscosity contrasts within the TBL are much smaller than values for the icy crust 
due to the viscosity-reducing effects of pressure and the shallower thermal gradient due to 
the melting curve. 



As for the depth-dependent parameters a(P, T) and Cp(T), we use results obtained 



for Ice VII. Fei et al. (1993) have produced the following expression for ct(P, T) , in units 
of K~ l : 



K' 



OT D \-0.9 



a(P, T) = (3.9 x 1(T 4 + 1.5 x KT b T)(l + -=^-P) 



(14) 



We again use the results from Fei et al. (1993) to obtain Cp in units of J kg K 



-i v-l. 



C P (P, T) = 4027.22 + 0.168T - 8.011 x 10 7 



J^2 



(15) 



Finally, we test throughgoing convection at phase transitions throughout the ice mantle. 
A sufficiently endothermic phase transition (i.e. the Clapeyron slope, 7, is sufficiently 
negative) can partition convection by introducing temperature-driven topography to the 
phase boundary. We evaluate the dimensionless "phase buoyancy parameter" (Pj) for all 



endothermic transitions, given by (Christensen and Yuen 1985): 
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p __Rb Ap 

Pb = 7-5- = i^rir (16) 

Ka p z gha 



Where Kb is the boundary Rayleigh number which uses a phase-driven instead of 
thermal-driven density contrast, 7 is the Clapeyron slope, Ap is the density contrast across 
the transition, p is the mean density, and h is the full thickness of the unpartitioned 
convective layer. Our values for Ap are calculated from EOSs for the relevant ice phases, 
while 7 was obtained directly from the literature or calculated with published triple point 
data. See Table H] for a listing of relevant references. 



4. Results 

We obtain mass-radius relationships for water planets from our hydrostatic models. 



Following the work of Valencia et al. (2007b), we fit our results to a function of total mass 
of the form: 

^<^> B < 17 > 

' Earth lvl Earth 

With values for A and B given for a range of water compositions in Table [2j The 
silicate - metal ratio in each case is 2.09. We find that Ice X lies at the bottom of the H 2 
shell of the 25% and 50% H 2 planets except for the 2M.Earthi 25% H 2 case (see Table [3]). 
In the cases of Earth-like water mass fraction, we find either Ice Ih or liquid water at the 
H 2 - silicate boundary, depending on the near-surface temperature profile (see below). 

For the planets with ice phase transitions (25% and 50% H 2 planets), our analysis of 
convective breakthrough revealed no phase transitions that should partition convection. As 
an example, the relevant parameters and P& values for the 5MEarth, 50% H 2 exoplanet 
(h ~ 4400 km) are listed in Table 111 Note that our predictions of convective breakthrough at 
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the Ice Ih - III and Ice II - V boundaries diverge from the corresponding results for the icy 



satellites (McKinnon 1998). This is due to the much greater thickness h of the convective 
cell in our exoplanet. For a much smaller exoplanet closer to the Galilean satellite instead 
of super-Earth regime, given an ice mantle thickness on the order of a few hundred km, 
these two transitions may hinder convection. 

Our models of the frozen surface case reveal that four qualitatively distinct dynamical 
Regimes exist for the near-surface region. For a given rheology, the Regimes are determined 
by the surface temperature Tg, heat flux q$, and, to a lesser extent, gravity g. The 
parameter space in which one finds each of the four layering regimes is shown for the 
g = 15 m/s 2 and g = 10 m/s 2 cases in Figure [3j A liquid ocean is present for Regimes I 
and II. In Regime I, the Ice Ih layer overlies a liquid water ocean and is not of sufficient 
thickness to convect. In Regime II, the temperature profile still intersects the melting 
curve within the Ice Ih region, thus maintaining the presence of an ocean. However, the 
Ice Ih region in this case is convectively unstable, and some solid-state convection is found 
at the base of this region. The boundary between Regime I and Regime II is calculated 
by setting the Rayleigh number of the strictly Ice Ih region above the ocean layer to the 
critical Rayleigh number (Ra crit ): if the Rayleigh number is above Ra cr a, the planet is in 
Regime II. The kink in this curve at m 225 K is due to a transition between the usage of 
stagnant-lid formulation for lower T$ and small viscosity contrast formulation for higher 
values. At lower T$ values, increasing the temperature lowers the viscosity and expands the 
parameter space where convection occurs. However, as T$ continues to increase, eventually 
the thinning of the crust inhibits convection, hence the characteristic "humped" shape of 
the Regime I - II boundary. 

In Regime III, the "turn off point" in the temperature profile- the point at which 
the switch is made from conductive to convective heat transfer- is below the minimum 
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temperature of the melting curve (251 K). Therefore, we find no liquid ocean. Instead, 
given the penetration of the Ice Ih - III and Ice Ih - II interfaces, we find a convection 
cell immediately below the crust that extends down to the silicate - ice boundary. Finally, 
in Regime IV, a thick, convectively stable crust extends below the region occupied by Ice 
Ih and meets the melting curve before the onset of convection. The temperature profile, 
when it reaches the melting temperature, is bounded by the melting curve. This must be 
the case: the temperature cannot be higher than the melting curve since the liquid would 
rapidly carry away heat, refreezing the material. The temperature profile must therefore 
follow the melting curve until this transition layer becomes convectively unstable, at which 
point it dives to greater depths along an adiabat. The dynamical properties of this region 
bound to the melt curve are considered further below under Discussion. 

The curve separating Regime III from Regimes II and IV is the contour for a constant 
"turn off point" temperature of 251 K. Meanwhile, the line separating Regimes I and 
IV corresponds to the scenario where the conductive portion of the temperature profile 
intersects the melt curve at the Liquid - Ice Ih - Ice III triple point: T = 251 K and 
P = 210 MPa. If the intersection with the melt curve occurs at lower pressure (above the 
Regime I - IV boundary in Figure [3]), the temperature profile intersects the melting curve 
within the Ice Ih domain and therefore falls into Regime I. 



Our formulation for Newtonian viscosity (Equation |4j) is subject to uncertainty due to 
the unknown characteristic grain size in the Ice Ih crust of exoplanets. Compared to the 
pa 0.2 mm we used in our analysis, Galileo observations of Europa has produced estimates of 
0.1mm for surficial grain size and estimates of several millimeters in the underlying regions 



(Pappalardo et al 1998). Assuming that volume diffusion continues to be the dominant 
creep mechanism, increasing the grain size to 0.5 mm increases the constant r] Q in Equation 
4] by a factor of six (rj oc d 2 ) and lowers the required qs value for all Regime boundaries 
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by about a factor of two, although their T$ coordinate remains essentially identical (see 
Figure [3]). The uncertainty in grain size, and hence Ice Ih viscosity, is the single greatest 
uncertainty in our analysis of the ice crust, and care should be taken to consider a range of 
possible grain sizes when attempting to characterize the surface of a specific water planet. 

A further source of uncertainty lies in the possibility of the presence of solutes or other 
volatiles in significant quantities within the ocean layer, thereby depressing the freezing 
point of the water-solute mixture, resulting in thinner and colder crusts, possibly hindering 
convection. We test the case of a water-ammonia mixture by adopting the phase diagram of 



the H2O - NH3 system as proposed by Hogenboom (1997). Assuming NH3 concentration 



of 10%, we find a modest correction to the Regime I - II boundary curve (see Figure [3]). 

As pointed out by the referee, we ignore the transition regime between the low viscosity 
contrast and stagnant-lid convection regimes. Given that our calculations yield Rayleigh 
numbers on the order of 10 5 to 10 6 for the Ice Ih shell in Regime II and > 10 6 for Regimes 
III and IV (equal to Ra of the full ice mantle convection in these latter cases), the surface 
temperature range over which the transitional regime is valid is between T5 rj 190 K and 
235 K. As the plausible deviation from our interpolation inside this narrow range is small 
compared to the magnitude of the other uncertainties that exist, we accept our current 
results as sufficiently accurate while noting that fuller consideration of the transitional 
regime may be worthwhile in a future work. 

Although it is not plotted in Figure [3j we imagine a separate regime where T$ is above 
freezing, resulting in a liquid ocean at the surface. The temperature profile in such an ocean 
is adiabatic until it freezes into higher phase ices at a depth on the order of 10s to a few 
100s of km. 

We apply our convection model for the ice mantle underlying a liquid ocean to the 25% 
and 50% H 2 exoplanets with T$ = 250 K. The near surface structure is first determined 
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in each case, yielding the result the crust in all these cases fall into Regime I. As an 
example, in the 5 M Earth, 50% H 2 case we find that an ocean of 62 km thickness underlies 
a brittle crust only 1.2 km thick. The ocean bottom temperature is Tt = 283 K; this 
figure is virtually constant for all six cases to within two degrees. The full solution to the 
temperature profile in the ice mantle for selected exoplanets is shown in Figure |4} We 
find that due to the smaller adiabatic temperature change and the correspondingly lower 
temperature at the top of the BBL, ATbbl is significantly larger for planets with the 
smallest H 2 shells. 

We can estimate the full ice mantle Rayleigh number of each planet. As described 
further in the Discussion section, finding a characteristic viscosity for the ice mantle is 
elusive due to its complex depth-dependence, reaching a viscosity maximum in mid-layer. 
We employ a simple approximation by taking the logarithmic average of the the viscosity 
throughout the ice mantle. This estimate results in Rayleigh numbers for the ice mantle 
between 10 5 and 10 9 , with the smallest Rayleigh numbers corresponding to the for the least 
massive planets with the thinnest water shell. 

The overturn time of the ice mantle (i.e. time needed for material from the silicate - ice 



boundary to reach the top of the ice mantle) can now be calculated. Turcotte and Schubert 



| (2002) have found by theoretical considerations that, given a cell aspect ratio of unity, 



the maximum velocity of convective material should depend on the full mantle Rayleigh 
number Ra m : 



u = 0.271 jRa 2 ^ (18) 

The overturn time scale is then simply calculated by dividing the thickness of the ice 
mantle by Uq. The results of this calculation for a range of planetary masses and water mass 
fractions, as well as their uncertainties, are shown in Table [5] A word about the associated 
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errors is found in Discussion. 



5. Discussion 

If a water planet has a frozen surface, the Regime of the crust has direct impact on 
observable features. In Regimes I and II, the presence of a liquid ocean layer is expected 
to bias the set of outgassed compounds towards species that are water-transportable, 
since, although solid-state convection through the mantle of higher ice phases can dredge 
up material expelled from the silicate - metal core, predominantly materials that can be 
dissolved or suspended in liquid water are able to traverse the ocean layer. A more detailed, 
quantitative analysis of this selection process is left for future work. 

In Regimes II, III, and IV, the presence of solid-state convection immediately below 
the brittle crust may lead to the existence of ice tectonics. This provides a means by which 
the chemistry of the convection cells can continuously gain a window to the surface, altering 



the atmospheric observables of the exoplanet (Valencia et al. 2007c). A simple test for 



the presence of ice tectonics is to compare the tensile strength of Ice Ih as found on the 



surfaces of icy satellites IMPa, Nimmo and Schenk 2006) and the expected tensile 



stress imposed on the crust by the convection cell. The approximate stress imposed by 



underlying convection is (McKinnon 1998): 

gaphAT a 



?>RaJ rit 



(19) 



Where again h and AT a are the height and driving temperature contrast of the 
convective region, respectively. The tensile stress on the brittle crust may then be amplified 
if the thickness of the crust D is smaller than the horizontal length over which a convection 
cell is in contact with the crust. Assuming a convection cell aspect ratio of unity, the stress 
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on the crust is therefore r = <j conv h/D. 



Performing this analysis for a crust in Regime II, we find that even in the case of small 
D (3 km) and high AT a (40 K), r does not rise beyond 220 kPa for g between 10 and 
20m/s 2 . Therefore, the crust in Regime II is probably not directly broken apart by stresses 
due to solid-state convection alone. However, resurfacing may persist in Regime II and 
even Regime I crusts due to a range of other phenomena that have been observed on the 
large icy satellites of the solar system. Such surface processes include diapirism (Europa- 



Pappalardo et al 1998; Triton- Schenk and Jackson 1993), cryovolcanism (Titan- Lopes et 



al 2007), and planetary-scale tectonics (Ganymede- Head et al 2002, Europa- Showman 



and Malhotra 1999), although the last of these maybe require exogenic sources of stress 



such as tidal forces in addition to internal processes. The presence of solid-state convection 



is found to facilitate these forms of resurfacing (e.g. diapirism- Pappalardo et al 1998 



cryovolcanism- Mitri et al. 2008), although penetration of the icy crust by underlying fluids 



may still be possible in the absence of convection (Showman et al 2004). Alternatively. 



a Regime II crust may be completely devoid of surface activity as in the case of Callisto 



(Showman and Malhotra 1999). 



As considerations of these solar system bodies show, reliably constraining the outgassing 
into the atmospheres of Regimes I and II planets is not possible given only the dynamical 
and thermal structure of the ice shell, as presented in this paper. 

On the other hand, in Regimes III and IV, the stress on the lithosphere is greatly 
increased for two reasons. First, because the convection cells extend to the base of the 
H2O layer of the planet, a conv is much larger due to greater values for h and AT CT . Second, 
assuming convective cells with aspect ratios of near unity, the amplification from the h/D 
factor is much greater. We find that in the 25% and 50% water cases considered in this 
work, the crust is easily broken by convection-driven stresses (r w 10 3 MPa for a thick, 
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50 km crust), leading to continuous global resurfacing under these Regimes. This result 
is consistent with analogous analysis for super-Earths with rocky surfaces; as shown in 



Valencia and O'Connell (2009), the higher convective stresses in the mantles of large 



silicate planets likewise outcompetes the tensile strength of the brittle crust. 

Finally, in the liquid surface case, chemical species released into the atmosphere should 
again be selected for water-transportable types, but outgassing on the surface proceeds 
unhindered. In this respect, the properties of out gassed material may be similar to those of 
Regimes I and II, while the flux may be similar to that of Regimes III and IV. 

Several uncertainties exist for our model of convection in the ice mantle. The true 
viscosity of the Ice VII and Ice X is unconstrained, as ice viscosity has not been measured 
at pressures of greater than about 0.8 GPa. The viscosity deep in the ice mantle is 
determined by the balance between three processes: increasing pressure, which increases 
the activation energy; shrinking activation volume, which tempers the above effect of 
pressure; and increasing temperature. Interestingly, all resulting mantle thermal profiles 
show a viscosity maximum in mid-layer of the ice mantle with viscosity values between 
10 21 and 10 23 Pa ■ s. This is due to the initial dominance of the pressure effect over the 
adiabatic temperature increase owing to a relatively large activitation volume. However, 
as the vacancies of the material become increasingly compressed at depth, the viscosity 
decreases with increasing temperature. This effect leads to low BBL viscosities of as 
little as 10 14 Pa ■ s for the 10 M Earth, 50% H 2 planet. Such results may not reflect true 
values and improved results must await the availability of additional rheological data for 
high-pressure ices. Furthermore, more refined results call for new theory that characterizes 
convection with a viscosity maximum in mid-layer and provides for a means to determine a 
characteristic viscosity. 

Values for certain other parameters in the mantle Rayleigh number, such as g and 



- 22 - 



a are difficult to choose to be representative of the ice mantle. We find that varying 
these throughout their possible range generally leads to a change of a factor of a few 
in the resulting overturn time. Combined with the greater uncertainty resulting from 
the characterization of viscosity as described above, the results in Table [5] are presented 
as estimates to within around one order of magnitude instead of precise solutions. We 
note, however, that despite these uncertainties, the qualitative phenomenon of a mid-layer 
high viscosity region is relatively certain, as the effect is apparent even when varying the 



parameters in the viscosity law beyond the range of published values (Durham et al. 1997). 



A further point of ambiguity is whether the steep rise in temperature along the adiabat 
in the lower parts of the ice mantle will cause our thermal profile to intersect the melting 
curve again at depth. Recent studies of the melting curve of high pressure ice phases 



suggest that this may occur only for a planet with a much larger water layer. Schwager 



et al. (2004) found experimentally that at 100 GPa corresponding to a the melting point 
of the Ice X is 2400 K, while our ice mantle temperatures are confined to under 1200 K 
(Figure EJ. 



A word must be said about the behavior of material where the temperature is 
constrained to the melting curve. This situation arises at the top of the ice mantle when 
the crust falls into Regimes I, II, or IV. In Regimes I and II, a liquid ocean overlies the ice 
mantle, and the temperature profile runs along the melt curve until convective instability is 
reached in this boundary layer. 

Although the globally averaged temperature proceeds along the melt curve in these 
cases, the full 3-D scenario is more complicated. Qualitatively, we expect the upward 
heat flux in this transition zone to be carried by the percolation of liquid water in certain 
concentrated zones. We describe a feedback mechanism that leads to the creation of 
narrow zones of rising meltwater. Beginning with a laterally homogeneous transition layer, 
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if a packet of material should melt in a specific location, say due to the presence of an 
underlying hot plume, this packet of material immediately undergoes a decrease in density, 
which leads to a decrease in the pressure in the underlying column. A decrease in pressure 
leads to the melting of ice previously at equilibrium, further promoting the creation of 
meltwater in the column. More detailed considerations of this process should be performed 
in the future. 

6. Conclusions 

In this work we have modeled the interior dynamics of the H 2 layer of water planets. 
We first study the thermal and dynamical properties of an ice crust for a planet with 
sub-freezing T$. We find that the dynamical behavior of the crust falls into one of four 
Regimes (Figure [3]), primarily determined by the planetary parameters surface temperature 
and heat flux, with a weaker dependence on surface gravity. Planets with a crust in 
Regimes III and IV should show atmospheric signatures of close to the full range of possible 
compounds produced from materials released from the silicate - metal core. Ice tectonics 
under these Regimes is expected to bring about continual resurfacing. The atmosphere 
of Regimes I and II planets may or may not be enriched by outgassing, depending on 
the presence of resurfacing phenomena such as cryovolcanism and other sources of stress 
(e.g. tides). As evidenced by the examples of such bodies studied in our solar system, the 
surface geology of icy bodies in these Regimes is difficult to constrain without case-specific 
considerations of a wider range of resurfacing processes. 

In the case of Regimes I and II, the range of materials that is able to reach the surface 
is affected by the presence of a liquid ocean layer, which selects for materials that can exist 
in solution or suspension. A more detailed study of this process is important for work future 
work. In the case of a liquid surface layer, we expect a full flux of outgassed materials 
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biased by its transportability through the liquid ocean. 

Our model of the ice mantle composed of higher-phase ices shows that materials 
expelled from the silicate - metal core into the ice layer can traverse the ice mantle on the 
time scale of 0.1 to a hundred Myr for the cases examined with shorter overturn times 
corresponding to larger planets. New constraints on the rheology of higher-phase ices as 
well as theory to describe the behavior of convection with a mid-layer viscosity maximum 
are required to obtain more accurate results. The specific chemical composition and mass 
flux of such dredged and outgassed material should be investigated quantitatively in further 
studies, ultimately leading to specific predictions of the atmospheric detectables of water 
planets. 
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Fig. 1. — The H 2 region of the 5 M^arth^ 50% H 2 planet drawn to scale. This configuration 
estimates heat flux according to Equation [2j 
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Fig. 2. — The qualitative features of the thermal profile within the ice mantle. The temper- 
ature in the TBL is bound to the melting curve while that of the BBL is conductive. 
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Fig. 3. — The four Regimes of the crust as a function of Ts and q$ for the g = 15 m/ s 2 and 10 m/s 2 cases. 
Contours of constant crustal thickness for g — 15 m/s 2 are plotted as light dashed lines in the top diagram. 
See Section [4] for explanation of the boundaries of each regime. Depth scale in the phase diagrams are based 
on the Earth, 50% H2O planet. The cartoon cross-sections are not drawn to scale. In the liquid surface 
case, the cross-section would appear as in Regime I, only without the "Crust" section; its thermal profile 
would begin at an above-freezing temperature at the surface and be adiabatic throughout the liquid ocean. 



- 33 - 



1000 
900 
800 



g 700 

I 

a 600 

i 

H 500 



400 



1 1 1 1 1 1 1 

- 




1 1 

- 




10 M Earlh , 50% H 2 




- 


5 M Earih , 50% H,0 J 






^^-rTT~25% H,0 
/ 5 M^/^>J 

^-^Z "" 2 M Earth , 25% 1 1 2 








■ 



0.0 0.2 0.4 0.6 0.8 1.0 

Proportion of ice mantle from ocean interface 



Fig. 4. — Complete temperature profile of the ice mantle of four selected water planets 
plotted against normalized depth from the ocean - Ice VI contact. True mantle thicknesses 
are: 1910 km, 2460 km, 4500 km, and 5425 km for the four planets plotted in ascending order 
as seen in the figure (i.e. 1910 km corresponds to the 2 MEarth, 25% H%0 planet, etc). The 
locations of the bottom boundary layers have been staggered in order to show the magnitude 
of temperature change. 
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Component 


EOS 


Po (^) 


Kqt (GPa) 


K' 


K^ T (GPa- 1 ) 


Ref. 


Fe 


Vinet 


8300 


156.2 


6.08 




1 


MgSi0 3 


B-M 4 


4100 


247 


3.97 


-0.016 


2 


IceX 


B-M 3 


1239 


4.26 


7.75 




3 


Ice VII 


B-M 3 


1463 


23.7 


4.15 




4 


Ice Ih 


B-M 3 


917 


9.86 


6.6 




5. 6 



Table 1: Equations of state and elastic parameters used in calculation of static planetary 
dimensions. B-M 3 and B-M 4 refers to the 3rd and 4th order Birch-Murnaghan equations of 
state while po an d Kq T are the zero-pressure density and pressure derivatives of the isothermal 
bulk modulus, repectively. 

References: (1) Anderson et al 2001| (2) |Karki and Wentzcovitch" 2000 (3) Loubeyre et al 



1999 (4) Hemley et al 1987] (5) CRC Handbook 1969(@ (6) Strassle et al. 2005 



Constant 


0.02% H 2 


25% H 2 


50% H 2 


A 


0.994 


1.157 


1.255 


B 


0.266 


0.253 


0.251 



Table 2: Constants for the mass-radius scaling relation (Equation 17) 



H 2 Content 


2 MEarth 


5 MEarth 


10 M E arth 


25% H 2 


47.8 ±1.0 GPa 


115.3 ± 1.6 GPa 


236.6 ±2.6 GPa 


50% H 2 


86.2 ±2.0 GPa 


212.0 ±2.2 GPa 


436.9 ±3.6 GPa 



Table 3: Pressures at the ice - silicate boundary. Values greater than about 60 GPa indicate 
the prescence of Ice X. Errors are determined by using the extremes of each equation of state 
parameter that yields the lowest and highest pressure. 
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Transition 


Ap (kg/m 3 ) 


7 (MPaK^ 1 ) 


Pb 


Ref. for Ap 


Ref. for 7 


Ih - III 


218 


-0.47 - -0.24 


-0.01 - -5 x 10~ 3 


1 


2; 3 


II -V 


38 


-8.9 - -0.67 


-0.025 - -.002 


1,4 


2; 3,5 


II -VI 


48 


« -1.7 


« -5 x 10~ 3 


1,4 


6 


VII -X 


68.1 


-11 - +26.3 


-0.01 - +0.02 


7,8 


9; 10 



Table 4: The phase buoyancy parameter for all encountered endothermic phase transitions 
between the water ices. Throughgoing convection is unhindered if P^ > —0.3. Since the 



Clapeyron slope 7 is the most uncertain parameter in Pj, (see Equation 16), the range of 
values of Pf, is found by evaluation with a range of 7 values found in the literature. 



References: (1) ]Shaw |p86| (2) |McKinnon |p98| (3) |Bridgman |p!2| (4) |Leon et al. ][2002 
(5) |Mercury et al~||2001| (6) |Durham et al~||1997| (7) |Loubeyre et al. ||1999| (8) |Hemley et al 
19871 ( 9 ) |Song et al. ||2003| (10) |Goncharov eraT~|p^999 



H 2 Content 


2 M Earth 


5 MEarth 


I® M E arth 


25% H 2 


83 Myr 


16 Myr 


0.8 Myr 


50% H 2 


190 Myr 


12 Myr 


0.3 Myr 



Table 5: Order of magnitude estimates for the overturn time scale for the ice mantles. Num- 
bers represent the time for material to rise from ice - silicate boundary to the top of the ice 
mantle. The values for the 10 MEarth planets are not as well-constrained. See the Discussion 
section for further notes about uncertainties. 



